
# Replication File

# All models are ran in the first section of this file
# Plots are generated in the second section


# R Version 3.4.2

# Packages versions
# Amelia 1.7.4
# plm 1.5-5
# lmtest 0.9-35

library(stargazer)
library(plm)
library(lmtest)
library(Amelia)

# Load data
data <- read.csv('fiscal_manipulation_data.csv',header=T)

# Variable Names (See SI 1.1 for definitions and data sources)

# log_con_inv: 'Consumption to Investment (Log)',
# prop_elapsed: 'Electoral Cycle',
# max_pm_gov_leg: 'Dissolution Powers',
# change_GDP: 'GDP Growth',
# log_gdp_dol: 'GDP (Log)',
# inflation: 'Inflation',
# revenue: 'Revenue',
# gross_debt: Debt to GDP',
# vote_share: 'Vote Share',
# govfrac: 'Government Frac.',
# log_con_inv_t: 'Consumption to Investment plus Transfers'
# log_gen_exp: 'General Gov. Expenditures'
# log_cent_exp: 'Central Gov. Expenditures'
# structural_bal: 'Cyclically Adjusted Primary Balance'
# transparency: 'Fiscal Transparency'
# BBR: 'Balanced Budget Rule'
# DR: 'Debt Rule'
# legislative_type: 'Electoral System'
# frac: 'Legislative Frac.',
# single_party: 'Single Party Gov.'
# execrlc: 'Partisan Orientation of Exec.',
# union_density: 'Union Density'
# max_elf: 'Ethno-Linguistic Fractionalization'
# dem_age: 'Democratic Age'
# fisc_cent: 'Fiscal Decentralization'
# eu: 'European Union'
# openness: 'Trade Openness'


# Impute data

set.seed(1234)

data_imp <- amelia(x = data, m=10, ts= "Year",
                   cs="country",ords=c('execrlc','eu','single_party','dem_age',
                                       'BBR','DR','legislative_type'), polytime = 3)

# Warnings shows high collinearity with debt rules and balanced budget rules, as noted in the MS
# Removing one removes collinearity (and warning) and 
# doesn't change results substantiallly




############################################
############ MODELS  ######################
########## (plots below) ##################
############################################

# Make function that takes reg equation as argument, outputs results

imp_fun <- function(form=y~x, mod='within',eff='time',data){ # note dat is changed
  
  b_out1 <- NULL # Ests
  se_out1 <- NULL # SEs
  r_square1 <- NULL # R2
  adj_r_square1 <-NULL # Adj R2

  for(i in 1:data$m){
    imp_data <- data$imputations[[i]]
    # Cubic polynomial to account for trend
    imp_data$t <- imp_data$Year - min(as.numeric(paste(imp_data$Year)))
    imp_data$t2 <- imp_data$t^2 
    imp_data$t3 <- imp_data$t^3
    
    
    mod1 <- plm(form,model=mod,effect = eff,index=c('country','Year'),data=imp_data) 
    
    se1 <-coeftest(mod1, vcov=vcovHC(mod1,cluster="group"))
    b_out1 <- rbind(b_out1, mod1$coef)
    se_out1 <- rbind(se_out1, se1[,2])
    r_square1 <- rbind(r_square1, summary(mod1)$r.square[1])
    adj_r_square1 <- rbind(adj_r_square1, summary(mod1)$r.square[2])
  }
  
  
  # Combine results from imputations
  combined_results1 <- mi.meld(q=b_out1, se=se_out1)
  # Get z scores
  mi_z1 <- combined_results1$q.mi/combined_results1$se.mi
  # Get p-values
  mi_p1 <- 2*(1 - pnorm(abs(mi_z1)))
  # Combine the above (minus z-score)
  combined_final1 <- cbind(t(round(combined_results1$q.mi,3)), t(round(combined_results1$se.mi,3)), round(t(mi_p1),4))
  colnames(combined_final1) <- c("EST", "SE", "P")
  
  r2 <- sum(r_square1)/length(r_square1)
  adj_r2 <- sum(adj_r_square1)/length(adj_r_square1)
  
  
  return(list(combined_final1,r2,adj_r2, nrow(imp_data)))
  
}


# Table 1, MODELS 1 - 8
### Tables printed out in LaTeX form (for comparison)

mod1_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                      vote_share+govfrac, 
                    mod="within",eff="time",data=data_imp)
# Can also print any of the imp models for comparison
## i.e. run the following line:
mod1_imp

mod2_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                      vote_share+govfrac+t+t2+t3,
                    mod="pooling",data=data_imp)

#mod2_imp if comparing in the Console

mod3_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                      vote_share+govfrac+t+t2+t3, 
                    mod="within", eff="individual",data=data_imp)

mod4_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                      vote_share+govfrac, mod="within",eff="twoways",data=data_imp)

mod5_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP,
                    mod="within",eff="time",data=data_imp)

mod6_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+t+t2+t3, mod="pooling",data=data_imp)

mod7_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+t+t2+t3, mod="within",eff="individual",data=data_imp)

mod8_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP, mod="within",eff="twoways",data=data_imp)


# Make baseline, non imputed, models, replace coefficients 
# THIS IS JUST FOR STARGAZER TO PRINT OUT imputed data

# Add cubic polynomial to baseline data
data2 <- data
data2$t <- data2$Year - min(as.numeric(paste(data2$Year)))
data2$t2 <- data2$t^2 
data2$t3 <- data2$t^3

mod1_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP+log_gdp_dol+inflation+revenue+gross_debt+vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data2)
mod2_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP+log_gdp_dol+inflation+revenue+gross_debt+vote_share+govfrac+t+t2+t3,index=c("country","Year"), model="pooling", data=data2)
mod3_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP+log_gdp_dol+inflation+revenue+gross_debt+vote_share+govfrac+t+t2+t3, model="within",effect="individual",index=c("country","Year"), data=data2)
mod4_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP+log_gdp_dol+inflation+revenue+gross_debt+vote_share+govfrac, model="within",effect="twoways",index=c("country","Year"), data=data2)

mod5_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP, model="within",effect="time",index=c("country","Year"), data=data2)
mod6_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP+t+t2+t3, model="pooling",index=c("country","Year"), data=data2)
mod7_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP+t+t2+t3, model="within",effect="individual",index=c("country","Year"), data=data2)
mod8_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+change_GDP, model="within",effect="twoways",index=c("country","Year"), data=data2)


# Print models out for LaTeX
stargazer(mod1_base,mod2_base,mod3_base,mod4_base,
          mod5_base,mod6_base,mod7_base,mod8_base,
          coef=c(list(mod1_imp[[1]][,1]),list(mod2_imp[[1]][,1]),
                 list(mod3_imp[[1]][,1]),list(mod4_imp[[1]][,1]),
                 list(mod5_imp[[1]][,1]),list(mod6_imp[[1]][,1]),
                 list(mod7_imp[[1]][,1]),list(mod8_imp[[1]][,1])),
          se=c(list(mod1_imp[[1]][,2]),list(mod2_imp[[1]][,2]),
               list(mod3_imp[[1]][,2]),list(mod4_imp[[1]][,2]),
               list(mod5_imp[[1]][,2]),list(mod6_imp[[1]][,2]),
               list(mod7_imp[[1]][,2]),list(mod8_imp[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption to Investment (Log)',
                              'Electoral Cycle','Dissolution Powers',
                              'GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.',
                              't','t2','t3','Electoral Cycle times Dissolution Powers'))

#R2
paste(round(mod1_imp[[2]],3),round(mod2_imp[[2]],3),
      round(mod3_imp[[2]],3),round(mod4_imp[[2]],3),
      round(mod5_imp[[2]],3),round(mod6_imp[[2]],3),
      round(mod7_imp[[2]],3),round(mod8_imp[[2]],3),sep=' & ')
# Adj R2
paste(round(mod1_imp[[3]],3),round(mod2_imp[[3]],3),
      round(mod3_imp[[3]],3),round(mod4_imp[[3]],3),
      round(mod5_imp[[3]],3),round(mod6_imp[[3]],3),
      round(mod7_imp[[3]],3),round(mod8_imp[[3]],3),sep=' & ')




# Supplemental Appendix

## SI 3.1
#Include transfers with DV

# Baseline model for Stargazer
trans_base <- plm(log_con_inv_t~lag(log_con_inv_t)+prop_elapsed*max_pm_gov_leg+
                    change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                    vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data2)

# Model with imputed data
trans_imp <- imp_fun(log_con_inv_t~lag(log_con_inv_t)+prop_elapsed*max_pm_gov_leg+
                       change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                       vote_share+govfrac, 
                     mod="within",eff="time",data=data_imp)


#Exps as DV (use difference)
gen_exp_base <- plm(diff(log_gen_exp)~lag(diff(log_gen_exp))+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                      vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data2)

gen_exp_imp <- imp_fun(diff(log_gen_exp)~lag(diff(log_gen_exp))+prop_elapsed*max_pm_gov_leg+
                         change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                         vote_share+govfrac, 
                       mod="within",eff="time",data=data_imp)

gen_exp_imp <- imp_fun(diff(log_gen_exp)~lag(diff(log_gen_exp))+prop_elapsed*max_pm_gov_leg+
                         change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                         vote_share+govfrac, 
                       mod="within",eff="time",data=data_imp)


cent_exp_base <- plm(diff(log_cent_exp)~lag(diff(log_cent_exp))+prop_elapsed*max_pm_gov_leg+
                       change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                       vote_share+govfrac
                     ,index=c("country","Year"), model="within",effect="time", data=data2)

cent_exp_imp <- imp_fun(diff(log_cent_exp)~lag(diff(log_cent_exp))+prop_elapsed*max_pm_gov_leg+
                          change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                          vote_share+govfrac, 
                        mod="within",eff="time",data=data_imp)


bal_base <- plm(structural_bal~lag(structural_bal)+prop_elapsed*max_pm_gov_leg+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data2)

bal_imp <- imp_fun(structural_bal~lag(structural_bal)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac, 
                   mod="within",eff="time",data=data_imp)

bal_base2 <- plm(structural_bal~lag(structural_bal)+prop_elapsed*max_pm_gov_leg+
                   change_GDP+log_gdp_dol+inflation+
                   vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data2)

bal_imp2 <- imp_fun(structural_bal~lag(structural_bal)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+
                      vote_share+govfrac, 
                    mod="within",eff="time",data=data_imp)


stargazer(trans_base,gen_exp_base,cent_exp_base,bal_base,bal_base2,
          coef=c(list(trans_imp[[1]][,1]),list(gen_exp_imp[[1]][,1]),
                 list(cent_exp_imp[[1]][,1]),list(bal_imp[[1]][,1]),
                 list(bal_imp2[[1]][,1])),
          se=c(list(trans_imp[[1]][,2]),list(gen_exp_imp[[1]][,2]),
               list(cent_exp_imp[[1]][,2]),list(bal_imp[[1]][,2]),
               list(bal_imp2[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption + Transfers to Investment (Log)','Gen. Gov. Expenditures',
                              'Cent. Gov. Expenditures','Adjusted Budget Balance',
                              'Electoral Cycle','Dissolution Powers',
                              'GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.',
                              'Electoral Cycle times Dissolution Powers'))


# R2
paste(round(trans_imp[[2]],3),round(gen_exp_imp[[2]],3),
      round(cent_exp_imp[[2]],3),round(bal_imp[[2]],3),
      round(bal_imp2[[2]],3),sep=' & ')
# Adj R2
paste(round(trans_imp[[3]],3),round(gen_exp_imp[[3]],3),
      round(cent_exp_imp[[3]],3),round(bal_imp[[3]],3),
      round(bal_imp2[[3]],3),sep=' & ')

paste(rep(trans_imp[[4]],5),sep=' & ')







## SI 3.2: (Changes to IV)

## NoteL Two imputed models that don't use function
  # 1) Gets subset to before 2012 (when DP index stops)

  # 2) Gets remoes cases after change in DP

# Two shorter tables; then no interaction term

b_out1 <- NULL # Ests
se_out1 <- NULL # SEs
r_square1 <- NULL # R2
adj_r_square1 <-NULL # Adj R2
for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  dat$t <- dat$Year - min(as.numeric(paste(dat$Year)))
  dat$t2 <- dat$t^2
  dat$t3 <- dat$t^3
  
  dat <- subset(dat,Year<=2011)
  
  mod1 <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                vote_share+govfrac,index=c("country","Year"), 
              model="within",effect='time', data=dat)
  
  # Standard Errors
  se1 <-coeftest(mod1, vcov=vcovHC(mod1,cluster="group"))
  # Get estimates
  b_out1 <- rbind(b_out1, mod1$coef)
  se_out1 <- rbind(se_out1, se1[,2])
  # Get r2/adj r2
  r_square1 <- rbind(r_square1, summary(mod1)$r.square[1])
  adj_r_square1 <- rbind(adj_r_square1, summary(mod1)$r.square[2])

}

combined_results1 <- mi.meld(q=b_out1, se=se_out1)
# Get z scores
mi_z1 <- combined_results1$q.mi/combined_results1$se.mi
# Get p-values
mi_p1 <- 2*(1 - pnorm(abs(mi_z1)))
# Combine the above (minus z-score)

# here is 3.2a, model1
pre_imp <- cbind(t(round(combined_results1$q.mi,3)), t(round(combined_results1$se.mi,3)), round(t(mi_p1),4))
colnames(pre_imp) <- c("EST", "SE", "P")
pre_r2 <- sum(r_square1)/length(r_square1)
pre_adj_r2 <- sum(adj_r_square1)/length(adj_r_square1)
pre_n <- nrow(dat)

# Get baseline for stargazer to work
pre_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac,index=c("country","Year"), 
                model="within",effect='time', data=subset(data,Year<=2011))


# Now 3.2a model 2
b_out1 <- NULL # Ests
se_out1 <- NULL # SEs
r_square1 <- NULL # R2
adj_r_square1 <-NULL # Adj R2
for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  dat$t <- dat$Year - min(as.numeric(paste(dat$Year)))
  dat$t2 <- dat$t^2
  dat$t3 <- dat$t^3
  
  dat <- subset(dat, !(country=='Czech Republic'&Year>=2009)&
                  !(country=='Slovak Republic'&Year>=1999)&
                  !(country=='United Kingdom'&Year>=2011))
  
  mod1 <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                vote_share+govfrac,index=c("country","Year"), 
              model="within",effect='time', data=dat)
  
  # Standard Errors
  se1 <-coeftest(mod1, vcov=vcovHC(mod1,cluster="group"))
  # Get estimates
  b_out1 <- rbind(b_out1, mod1$coef)
  se_out1 <- rbind(se_out1, se1[,2])
  
  r_square1 <- rbind(r_square1, summary(mod1)$r.square[1])
  adj_r_square1 <- rbind(adj_r_square1, summary(mod1)$r.square[2])
  
}
combined_results1 <- mi.meld(q=b_out1, se=se_out1)
# Get z scores
mi_z1 <- combined_results1$q.mi/combined_results1$se.mi
# Get p-values
mi_p1 <- 2*(1 - pnorm(abs(mi_z1)))
# Combine the above (minus z-score)
# 3.2a model 2
change_imp <- cbind(t(round(combined_results1$q.mi,3)), t(round(combined_results1$se.mi,3)), round(t(mi_p1),4))
colnames(change_imp) <- c("EST", "SE", "P")
change_imp
change_r2 <- sum(r_square1)/length(r_square1)
change_adj_r2 <- sum(adj_r_square1)/length(adj_r_square1)
change_n <- nrow(dat)

# Non-imputed model for stargazer
change_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), 
                   model="within",effect='time', data=subset(data, !(country=='Czech Republic'&Year>=2009)&
                                                               !(country=='Slovak Republic'&Year>=1999)&
                                                               !(country=='United Kingdom'&Year>=2011)))

# 3.2a model 3: no interaction (Use function defined above)
no_int<- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed+max_pm_gov_leg+
                   change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                   vote_share+govfrac, 
                 mod="within",eff="time",data=data_imp)


# Make baseline models, replace coefficients
no_int_base<- plm(log_con_inv~lag(log_con_inv)+prop_elapsed+max_pm_gov_leg+change_GDP+log_gdp_dol+inflation+revenue+gross_debt+vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data2)


# Now print the model for LaTeX
stargazer(pre_base,change_base,no_int_base,
          coef=c(list(pre_imp[,1]),list(change_imp[,1]),
                 list(no_int[[1]][,1])),
          se=c(list(pre_imp[,2]),list(change_imp[,2]),
               list(no_int[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption to Investment (Log)',
                              'Electoral Cycle','Dissolution Powers',
                              'GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.',
                              'Electoral Cycle times Dissolution Powers'))


# R2
paste(round(pre_r2,3),round(change_r2,3),round(no_int[[2]],3),sep=' & ')

# Adj R2
paste(round(pre_adj_r2,3),round(pre_adj_r2,3),round(no_int[[3]],3),
      sep=' & ')

pre_n
change_n
no_int[[4]]


##########################################
############## SI 3.2B ##############

# transparency; fiscal rules; electoral systems

# Fiscal transparency

# model 1
transp_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+transparency,index=c("country","Year"), model="within",effect="time", data=data)

transp_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+govfrac+transparency, 
                      mod="within",eff="time",data=data_imp)

# Fiscal Rules

# balanced budget rule (Model 2)
BBR_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac+BBR,index=c("country","Year"), model="within",effect="time", data=data)


BBR_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+BBR, 
                   mod="within",eff="time",data=data_imp)

# debt rule (model 3)
DR_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                 change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                 vote_share+govfrac+DR,index=c("country","Year"), model="within",effect="time", data=data)

DR_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                    change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                    vote_share+govfrac+DR, 
                  mod="within",eff="time",data=data_imp)


# Legislative type (mod 4) (Note: 1=maj, 2=PR, 3=mixed)

leg_type_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                       change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                       vote_share+govfrac+as.factor(legislative_type),index=c("country","Year"), model="within",effect="time", data=data)

leg_type_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                          change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                          vote_share+govfrac+as.factor(legislative_type), 
                        mod="within",eff="time",data=data_imp)


joint1_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+transparency+BBR+DR+as.factor(legislative_type),index=c("country","Year"), model="within",effect="time", data=data)

joint1_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+govfrac+transparency+BBR+DR+as.factor(legislative_type), 
                      mod="within",eff="time",data=data_imp)


stargazer(transp_base,BBR_base,DR_base,leg_type_base,joint1_base,
          coef=c(list(transp_imp[[1]][,1]),list(BBR_imp[[1]][,1]),list(DR_imp[[1]][,1]),
                 list(leg_type_imp[[1]][,1]),list(joint1_imp[[1]][,1])),
          se=c(list(transp_imp[[1]][,2]),list(BBR_imp[[1]][,2]),list(DR_imp[[1]][,2]), 
               list(leg_type_imp[[1]][,2]),list(joint1_imp[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption to Investment (Log)',
                              'Electoral Cycle','Dissolution Powers',
                              'GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.',
                              'Fiscal Transparency', 'Balanced Budget Rules','Debt Rules',
                              'Proportional Rep.','Mixed Rep.','Electoral Cycle times Dissolution Powers'))


# R2
paste(round(transp_imp[[2]],3),round(BBR_imp[[2]],3),
      round(DR_imp[[2]],3),round(leg_type_imp[[2]],3),
      round(joint1_imp[[2]],3),sep=' & ')

# Adj R2
paste(round(transp_imp[[3]],3),round(BBR_imp[[3]],3),
      round(DR_imp[[3]],3),round(leg_type_imp[[3]],3),
      round(joint1_imp[[3]],3),sep=' & ')

pre_n
change_n
leg_type_imp[[4]]





############################################################
############################################################
############### SI 3.2 C PLACEBO test

# 2a) Placebo
# Fiscal transparency

# Model 1
transp_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*transparency+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)

transp_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*transparency+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+govfrac, 
                      mod="within",eff="time",data=data_imp)

# balanced budget rule (Model 2)
BBR_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*BBR+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)


BBR_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*BBR+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac, 
                   mod="within",eff="time",data=data_imp)

# debt rule (model 3)
DR_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*DR+
                 change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                 vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)

DR_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*DR+
                    change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                    vote_share+govfrac, 
                  mod="within",eff="time",data=data_imp)


# Legislative type (Model 4) (Note: 1=maj, 2=PR, 3=mixed)

leg_type_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*as.factor(legislative_type)+
                       change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                       vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)

leg_type_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*as.factor(legislative_type)+
                          change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                          vote_share+govfrac, 
                        mod="within",eff="time",data=data_imp)

stargazer(transp_base,BBR_base,DR_base,leg_type_base,
          coef=c(list(transp_imp[[1]][,1]),list(BBR_imp[[1]][,1]),
                 list(DR_imp[[1]][,1]),list(leg_type_imp[[1]][,1])),
          se=c(list(transp_imp[[1]][,2]),list(BBR_imp[[1]][,2]),
               list(DR_imp[[1]][,2]),list(leg_type_imp[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption to Investment (Log)',
                              'Electoral Cycle','Transparency','BBR','DR',
                              'Right Exec.','Undefined Exec.','GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.',
                              'Elapsed times Transparency','Elapsed times BBR','Elapsed times DR',
                              'Elapsed times Proportional Rep','Elapsed times Mixed Rep.'))

#R2
paste(round(transp_imp[[2]],3),round(BBR_imp[[2]],3),
      round(DR_imp[[2]],3),round(leg_type_imp[[2]],3),sep=' & ')
# Adj R2
paste(round(transp_imp[[3]],3),round(BBR_imp[[3]],3),
      round(DR_imp[[3]],3),round(leg_type_imp[[3]],3),sep=' & ')


###################################################
################# SI 3.2D: Triple Interactions
###################################################

# Model 1
transp_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*transparency+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)

transp_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*transparency+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+govfrac, 
                      mod="within",eff="time",data=data_imp)

# Fiscal Rules

# balanced budget rule (Model 2)
BBR_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*BBR+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)


BBR_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*BBR+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac, 
                   mod="within",eff="time",data=data_imp)

# debt rule (Model 3)
DR_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*DR+
                 change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                 vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)

DR_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*DR+
                    change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                    vote_share+govfrac, 
                  mod="within",eff="time",data=data_imp)

# Legislative type (Model 4) (Note: 1=maj, 2=PR, 3=mixed)

leg_type_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*as.factor(legislative_type)+
                       change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                       vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=data)

leg_type_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*as.factor(legislative_type)+
                          change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                          vote_share+govfrac, 
                        mod="within",eff="time",data=data_imp)


stargazer(transp_base,BBR_base,DR_base,leg_type_base,
          coef=c(list(transp_imp[[1]][,1]),list(BBR_imp[[1]][,1]),
                 list(DR_imp[[1]][,1]),list(leg_type_imp[[1]][,1])),
          se=c(list(transp_imp[[1]][,2]),list(BBR_imp[[1]][,2]),
               list(DR_imp[[1]][,2]),list(leg_type_imp[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption to Investment (Log)',
                              'Electoral Cycle','Dissolution Powers','Transparency','BBR','DR',
                              'Proportional Rep.','Mixed Rep.','GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.', 
                              'Elapsed times Dissolution Powers', 'Elapsed times Transparency','Dissolution Powers times Transparency','Elapsed times Dissolution Powers times Transparency',
                              'Elapsed times BBR','Dissolution Powers times BBR','Elapsed times Dissolution Powers times BBR',
                              'Elapsed times DR','Dissolution Powers times DR','Elapsed times Dissolution Powers times DR',
                              'Elapsed times Proportional Rep.','Elapsed times Mixed Rep.','Dissolution Powers times Proportional Rep.','Dissolution Powers times Mixed Rep.',
                              'Elapsed times Dissolution Powers times Proportional Rep.','Elapsed times Dissolution Powers times Mixed Rep.'))

#R2
paste(round(transp_imp[[2]],3),round(BBR_imp[[2]],3),
      round(DR_imp[[2]],3),round(leg_type_imp[[2]],3),sep=' & ')
# Adj R2
paste(round(transp_imp[[3]],3),round(BBR_imp[[3]],3),
      round(DR_imp[[3]],3),round(leg_type_imp[[3]],3),sep=' & ')




###################################################
################# SI 3.3: Socio-Political Confounders
###################################################



# Model 1: Legislative frac
leg_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac+frac,index=c("country","Year"), model="within",effect="time", data=data2)

leg_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+frac, 
                   mod="within",eff="time",data=data_imp)

# Model 2: Single party

single_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+single_party,index=c("country","Year"), model="within",effect="time", data=data)


single_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+single_party, 
                      mod="within",eff="time",data=data_imp)


# Model 3: Union Density
union_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                    change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                    vote_share+govfrac+union_density,index=c("country","Year"), model="within",effect="time", data=data2)

union_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                       change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                       vote_share+govfrac+union_density, 
                     mod="within",eff="time",data=data_imp)

# Model 4: Ethno-Linguistic Fractionalization
elf_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac+max_elf,index=c("country","Year"), model="within",effect="time", data=data2)

elf_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+max_elf, 
                   mod="within",eff="time",data=data_imp)


# Model 5: Executive left-right orientation
exec_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                   change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                   vote_share+govfrac+as.factor(execrlc),index=c("country","Year"), model="within",effect="time", data=data2)

exec_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                      vote_share+govfrac+as.factor(execrlc), 
                    mod="within",eff="time",data=data_imp)

# Model 6: Full model (all above potential confounders)
joint2_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+frac+single_party+union_density+max_elf+as.factor(execrlc),index=c("country","Year"), model="within",effect="time", data=data)

joint2_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+govfrac+frac+single_party+union_density+max_elf+as.factor(execrlc), 
                      mod="within",eff="time",data=data_imp)




stargazer(leg_base,single_base,union_base,elf_base,exec_base,joint2_base,
          coef=c(list(leg_imp[[1]][,1]),list(single_imp[[1]][,1]),
                 list(union_imp[[1]][,1]),list(elf_imp[[1]][,1]),
                 list(exec_imp[[1]][,1]),list(joint2_imp[[1]][,1])),
          se=c(list(leg_imp[[1]][,2]),list(single_imp[[1]][,2]),
               list(union_imp[[1]][,2]),list(elf_imp[[1]][,2]),
               list(exec_imp[[1]][,2]),list(joint2_imp[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption to Investment (Log)',
                              'Electoral Cycle','Dissolution Powers',
                              'GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.',
                              'Leglislative Frac.','Single Party Gov.',
                              'Union Density','Ethno-Linguistic Frac.',
                              'Center Exec.','Right Exec.','Undefined Exec.',
                              'Electoral Cycle times Dissolution Powers'))


paste(round(leg_imp[[2]],3),round(single_imp[[2]],3),
      round(union_imp[[2]],3),round(elf_imp[[2]],3),
      round(exec_imp[[2]],3),round(joint2_imp[[2]],3), sep=' & ')
# Adj R2
paste(round(leg_imp[[3]],3),round(single_imp[[3]],3),
      round(union_imp[[3]],3),round(elf_imp[[3]],3),
      round(exec_imp[[3]],3),round(joint2_imp[[3]],3),sep=' & ')

elf_imp[[4]]

#######################################################################
################# SI 3.4: Macro Political and Economic Confounders
#######################################################################

# Model 1: Democratic age
age_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                  change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                  vote_share+govfrac+dem_age,index=c("country","Year"), model="within",effect="time", data=data2)


age_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+dem_age, 
                   mod="within",eff="time",data=data_imp)


# Model 2: European Union
eu_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                 change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                 vote_share+govfrac+eu,index=c("country","Year"), model="within",effect="time", data=data2)

eu_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                    change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                    vote_share+govfrac+eu, 
                  mod="within",eff="time",data=data_imp)

# Model 3: Fiscl decentralization
decent_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+fisc_cent,index=c("country","Year"), model="within",effect="time", data=data2)

decent_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+govfrac+fisc_cent, 
                      mod="within",eff="time",data=data_imp)

# Model 4: Trade Openness

open_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                   change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                   vote_share+govfrac+openness,index=c("country","Year"), model="within",effect="time", data=data2)

open_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                      change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                      vote_share+govfrac+openness, 
                    mod="within",eff="time",data=data_imp)

# Model 5: All above potential confounders
joint3_base <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac+dem_age+fisc_cent+eu+openness,index=c("country","Year"), model="within",effect="time", data=data2)

joint3_imp <- imp_fun(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                        change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                        vote_share+govfrac+dem_age+eu+fisc_cent+openness, 
                      mod="within",eff="time",data=data_imp)





stargazer(age_base,eu_base,decent_base,open_base,joint3_base,
          coef=c(list(age_imp[[1]][,1]),list(eu_imp[[1]][,1]),
                 list(decent_imp[[1]][,1]),
                 list(open_imp[[1]][,1]),list(joint3_imp[[1]][,1])),
          se=c(list(age_imp[[1]][,2]),list(eu_imp[[1]][,2]),
               list(decent_imp[[1]][,2]),
               list(open_imp[[1]][,2]),list(joint3_imp[[1]][,2])),
          omit.stat = c('f'),
          covariate.labels =c('Consumption to Investment (Log)',
                              'Electoral Cycle','Dissolution Powers',
                              'GDP Growth','GDP (Log)','Inflation',
                              'Revenue','Debt to GDP','Vote Share','Government Frac.',
                              'Democratic Age','EU','Fiscal Decentralization','Trade Openness',
                              'Electoral Cycle times Dissolution Powers'))


# R2
paste(round(age_imp[[2]],3),round(eu_imp[[2]],3),
      round(decent_imp[[2]],3),
      round(open_imp[[2]],3),round(joint3_imp[[2]],3),sep=' & ')
# Adj R2
paste(round(age_imp[[3]],3),round(eu_imp[[3]],3),
      round(decent_imp[[3]],3),
      round(open_imp[[3]],3),round(joint3_imp[[3]],3),sep=' & ')


decent_imp[[4]]






############################################
#################################
######## Plots
############################################

# Figure 1

dy_de_coef <- NULL
dy_dp_coef <- NULL
dy_de_se <- NULL
dy_dp_se <- NULL

for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  #  dat <- pdata.frame(dat, index=c("country", "Year"))
  
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share,index=c("country","Year"), 
                   model="within",effect='time', data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  dy_de_coef <- rbind(dy_de_coef, ests["prop_elapsed"] + ests["prop_elapsed:max_pm_gov_leg"]*DP)
  dy_dp_coef <- rbind(dy_dp_coef, ests["max_pm_gov_leg"] + ests["prop_elapsed:max_pm_gov_leg"]*E)
  
  dy_de_se <- rbind(dy_de_se, sqrt(cov_mat["prop_elapsed","prop_elapsed"] +DP^2*cov_mat["prop_elapsed:max_pm_gov_leg","prop_elapsed:max_pm_gov_leg"] + 2*DP*cov_mat["prop_elapsed","prop_elapsed:max_pm_gov_leg"]))
  dy_dp_se <- rbind(dy_dp_se, sqrt(cov_mat["max_pm_gov_leg","max_pm_gov_leg"] +E^2*cov_mat["prop_elapsed:max_pm_gov_leg","prop_elapsed:max_pm_gov_leg"] + 2*E*cov_mat["max_pm_gov_leg","prop_elapsed:max_pm_gov_leg"]))
}

combined_1 <- mi.meld(q=dy_de_coef, se=dy_de_se)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_dp_coef, se=dy_dp_se)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi

# Figure 1a

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="Marginal Effect of \n Proportion of Term Elapsed", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)


# Figure 1b
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main="Marginal Effect of \n Dissolution Powers", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)



################## SI PLOTS ####################################

# Write function to calculate ME of triple interaction
calc_dx <- function(x,z,w,est,Z,W){
  xz = paste(x,z,sep=':')
  xw = paste(x,w,sep=':')
  xzw = paste(x,z,w ,sep=':')
  return(
    est[x] + est[xz]*Z  + est[xw]*W + est[xzw]*Z*W)
}

calc_dz <- function(x,z,w,est,X,W){
  xz = paste(x,z,sep=':')
  zw = paste(z,w,sep=':')
  xzw = paste(x,z,w ,sep=':')
  return(
    est[z] + est[xz]*X  + est[zw]*W + est[xzw]*X*W)
}


calc_dx_se <- function(x,z,w,cm,Z,W){
  xz <- paste(x,z,sep=':')
  xw <- paste(x,w,sep=':')
  xzw <- paste(x,z,w ,sep=':')
  out <- sqrt(cm[x,x] + (Z^2)*cm[xz,xz] +(W^2)*cm[xw,xw]+(Z^2)*(W^2)*cm[xzw,xzw] +
                2*Z*cm[x,xz] + 2*W*cm[x,xw] + 2*Z*W*cm[x,xzw] + 2*Z*W*cm[xz,xw]+
                2*W*(Z^2)*cm[xz,xzw] + 2*Z*(W^2)*cm[xw,xzw]) 
  return(out
  )
}


calc_dz_se <- function(x,z,w,cm,X,W){
  xz = paste(x,z,sep=':')
  zw = paste(z,w,sep=':')
  xzw = paste(x,z,w ,sep=':')
  return(
    sqrt(cm[z,z] + (X^2)*cm[xz,xz] +(W^2)*cm[zw,zw]+(X^2)*(W^2)*cm[xzw,xzw] +
           2*X*cm[z,xz] + 2*W*cm[z,zw] + 2*X*W*cm[z,xzw] + 2*X*W*cm[xz,zw]+
           2*W*(X^2)*cm[xz,xzw] + 2*X*(W^2)*cm[zw,xzw]))
}

########## Fig. 3.2d1: Transparency ##########

dy_de_coef1 <- NULL
dy_de_coef2 <- NULL

dy_de_se1 <- NULL
dy_de_se2 <- NULL


for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  #  dat <- pdata.frame(dat, index=c("country", "Year"))
  
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*transparency+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  Transp <- c(0,50,100)
  
  dy_de_coef1 <- rbind(dy_de_coef1, calc_dx(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='transparency',est=ests,Z=DP, W = 50))
  
  dy_de_se1 <- rbind(dy_de_se1, 
                     calc_dx_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='transparency',cm=cov_mat,Z=DP, W = 50))
  
  dy_de_coef2 <- rbind(dy_de_coef2, calc_dz(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='transparency',est=ests,X=E, W = 50))
  
  dy_de_se2 <- rbind(dy_de_se2, 
                     calc_dz_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='transparency',cm=cov_mat,X=E, W = 50))
  
  
}

combined_1 <- mi.meld(q=dy_de_coef1, se=dy_de_se1)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_de_coef2, se=dy_de_se2)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi


# Fig 3.2d.1(a)

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="Transparency = 50", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)

# Fig 3.2d.1(b)
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main=" Transparency = 50", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)


########## Fig. 3.2d3: Balanced Budget Rules (set to zero)


dy_de_coef1 <- NULL
dy_de_coef2 <- NULL

dy_de_se1 <- NULL
dy_de_se2 <- NULL


for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  #  dat <- pdata.frame(dat, index=c("country", "Year"))
  
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*BBR+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  Transp <- c(0,50,100)
  
  dy_de_coef1 <- rbind(dy_de_coef1, calc_dx(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='BBR',est=ests,Z=DP, W = 0))
  
  dy_de_se1 <- rbind(dy_de_se1, 
                     calc_dx_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='BBR',cm=cov_mat,Z=DP, W = 0))
  
  dy_de_coef2 <- rbind(dy_de_coef2, calc_dz(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='BBR',est=ests,X=E, W = 0))
  
  dy_de_se2 <- rbind(dy_de_se2, 
                     calc_dz_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='BBR',cm=cov_mat,X=E, W = 0))
  
}

combined_1 <- mi.meld(q=dy_de_coef1, se=dy_de_se1)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_de_coef2, se=dy_de_se2)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi


# Figure 3.2d.2(a)

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="BBR = 0", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)


# Figure 3.2d.2(b)
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
#hist(data2$BBR, xlab=NULL, ylab=NULL, col = "light blue"
#     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main="BBR = 0", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)


########## Fig. 3.2d3: Balanced Budget Rules (set to one)

dy_de_coef1 <- NULL
dy_de_coef2 <- NULL

dy_de_se1 <- NULL
dy_de_se2 <- NULL


for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  #  dat <- pdata.frame(dat, index=c("country", "Year"))
  
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*BBR+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  Transp <- c(0,50,100)
  
  dy_de_coef1 <- rbind(dy_de_coef1, calc_dx(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='BBR',est=ests,Z=DP, W = 1))
  
  dy_de_se1 <- rbind(dy_de_se1, 
                     calc_dx_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='BBR',cm=cov_mat,Z=DP, W = 1))
  
  dy_de_coef2 <- rbind(dy_de_coef2, calc_dz(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='BBR',est=ests,X=E, W = 1))
  
  dy_de_se2 <- rbind(dy_de_se2, 
                     calc_dz_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='BBR',cm=cov_mat,X=E, W = 1))
  
}

combined_1 <- mi.meld(q=dy_de_coef1, se=dy_de_se1)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_de_coef2, se=dy_de_se2)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi


# Figure 3.2d.2(a)
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="BBR = 1", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)

# Figure 3.2d.2(b)
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
#hist(data2$BBR, xlab=NULL, ylab=NULL, col = "light blue"
#     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main="BBR = 1", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)


########## Fig. 3.2d4: Balanced Budget Rules (set to zero)

dy_de_coef1 <- NULL
dy_de_coef2 <- NULL

dy_de_se1 <- NULL
dy_de_se2 <- NULL


for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  #  dat <- pdata.frame(dat, index=c("country", "Year"))
  
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*DR+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  
  dy_de_coef1 <- rbind(dy_de_coef1, calc_dx(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='DR',est=ests,Z=DP, W = 0))
  
  dy_de_se1 <- rbind(dy_de_se1, 
                     calc_dx_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='DR',cm=cov_mat,Z=DP, W = 0))
  
  dy_de_coef2 <- rbind(dy_de_coef2, calc_dz(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='DR',est=ests,X=E, W = 0))
  
  dy_de_se2 <- rbind(dy_de_se2, 
                     calc_dz_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='DR',cm=cov_mat,X=E, W = 0))
  
}

combined_1 <- mi.meld(q=dy_de_coef1, se=dy_de_se1)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_de_coef2, se=dy_de_se2)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi


### Fig. 3.2d.4(a)
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="DR = 0", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)


### Fig. 3.2d.4(b)

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
#hist(data2$BBR, xlab=NULL, ylab=NULL, col = "light blue"
#     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main="DR = 0", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)


########## Fig. 3.2d5: Balanced Budget Rules (set to one)

dy_de_coef1 <- NULL
dy_de_coef2 <- NULL

dy_de_se1 <- NULL
dy_de_se2 <- NULL


for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  #  dat <- pdata.frame(dat, index=c("country", "Year"))
  
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*DR+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  
  dy_de_coef1 <- rbind(dy_de_coef1, calc_dx(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='DR',est=ests,Z=DP, W = 1))
  
  dy_de_se1 <- rbind(dy_de_se1, 
                     calc_dx_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='DR',cm=cov_mat,Z=DP, W = 1))
  
  dy_de_coef2 <- rbind(dy_de_coef2, calc_dz(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='DR',est=ests,X=E, W = 1))
  
  dy_de_se2 <- rbind(dy_de_se2, 
                     calc_dz_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='DR',cm=cov_mat,X=E, W = 1))
  
}

combined_1 <- mi.meld(q=dy_de_coef1, se=dy_de_se1)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_de_coef2, se=dy_de_se2)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi


### Fig. 3.2d.5(a)

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="DR = 1", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)

### Fig. 3.2d.5(b)

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
#hist(data2$BBR, xlab=NULL, ylab=NULL, col = "light blue"
#     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main="DR = 1", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)



########## Fig. 3.2d6: Proportional Representation

dy_de_coef1 <- NULL
dy_de_coef2 <- NULL

dy_de_se1 <- NULL
dy_de_se2 <- NULL


for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  dat$pr <- ifelse(dat$legislative_type==2,1,0)
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*pr+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  
  dy_de_coef1 <- rbind(dy_de_coef1, calc_dx(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='pr',est=ests,Z=DP, W = 1))
  
  dy_de_se1 <- rbind(dy_de_se1, 
                     calc_dx_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='pr',cm=cov_mat,Z=DP, W = 1))
  
  dy_de_coef2 <- rbind(dy_de_coef2, calc_dz(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='pr',est=ests,X=E, W = 1))
  
  dy_de_se2 <- rbind(dy_de_se2, 
                     calc_dz_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='pr',cm=cov_mat,X=E, W = 1))
  
}

combined_1 <- mi.meld(q=dy_de_coef1, se=dy_de_se1)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_de_coef2, se=dy_de_se2)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi



### Fig. 3.2d6(a)
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="Proportional Rep. = 1", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)


### Fig. 3.2d6(b)
par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
#hist(data2$BBR, xlab=NULL, ylab=NULL, col = "light blue"
#     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main="Proportional Rep. = 1", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)


########## Fig. 3.2d7: Mixed Representation


dy_de_coef1 <- NULL
dy_de_coef2 <- NULL

dy_de_se1 <- NULL
dy_de_se2 <- NULL


for(i in 1:data_imp$m){
  dat <- data_imp$imputations[[i]]
  #  dat <- pdata.frame(dat, index=c("country", "Year"))
  dat$mixed <- ifelse(dat$legislative_type==3,1,0)
  mod1_pool <- plm(log_con_inv~lag(log_con_inv)+prop_elapsed*max_pm_gov_leg*mixed+
                     change_GDP+log_gdp_dol+inflation+revenue+gross_debt+
                     vote_share+govfrac,index=c("country","Year"), model="within",effect="time", data=dat)
  ests <- coef(mod1_pool)
  cov_mat <- vcovHC(mod1_pool,cluster="group") # For FE
  DP <- seq(0,10)
  E <- seq(0,10)/10
  
  dy_de_coef1 <- rbind(dy_de_coef1, calc_dx(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='mixed',est=ests,Z=DP, W = 1))
  
  dy_de_se1 <- rbind(dy_de_se1, 
                     calc_dx_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='mixed',cm=cov_mat,Z=DP, W = 1))
  
  dy_de_coef2 <- rbind(dy_de_coef2, calc_dz(x='prop_elapsed',
                                            z='max_pm_gov_leg',w='mixed',est=ests,X=E, W = 1))
  
  dy_de_se2 <- rbind(dy_de_se2, 
                     calc_dz_se(x='prop_elapsed',z='max_pm_gov_leg',
                                w='mixed',cm=cov_mat,X=E, W = 1))
  
}

combined_1 <- mi.meld(q=dy_de_coef1, se=dy_de_se1)
upper1 <- combined_1$q.mi + 1.96*combined_1$se.mi
lower1 <- combined_1$q.mi - 1.96*combined_1$se.mi

combined_2 <- mi.meld(q=dy_de_coef2, se=dy_de_se2)
upper2 <- combined_2$q.mi + 1.96*combined_2$se.mi
lower2 <- combined_2$q.mi - 1.96*combined_2$se.mi


### Fig. 3.2d7(a)

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$max_pm_gov_leg, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$max_pm_gov_leg, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(DP), max(DP)), c(min(lower1), max(upper1)), xlab="Dissolution Powers",
     ylab="", main="Mixed Rep. = 1", type="n") # type="n" removes points
lines(DP, combined_1$q.mi, lwd=2)
lines(DP, lower1, lty = 3, lwd=2) # lwd makes lines thicker
lines(DP, upper1, lty = 3, lwd=2)
mtext("ME of Term Elapsed", side=2, line=2)
abline(h=0)


### Fig. 3.2d7(b)

par(mar=c(4,4.1,4.1,3.2)) # more space on right, less on left

# histogram
hist(data2$prop_elapsed, xlab=NULL, ylab=NULL, col = "light gray"
     ,lty=0, axes = F, main=NULL)
#hist(data2$BBR, xlab=NULL, ylab=NULL, col = "light blue"
#     ,lty=0, axes = F, main=NULL)
axis(4) # Put labels on right of plot
mtext("Observations", side=4, line=2)
rug(jitter(data2$prop_elapsed, amount=1)) #rugplot
par(new=TRUE) # To overlap ME plot on histogram

plot(c(min(E), max(E)), c(min(lower2), max(upper2)), xlab="Term Elapsed",
     ylab="", main="Mixed Rep. = 1", type="n") # type="n" removes points
lines(E, combined_2$q.mi, lwd=2)
lines(E, lower2, lty = 3, lwd=2) # lwd makes lines thicker
lines(E, upper2, lty = 3, lwd=2)
mtext("ME of Diss. Powers", side=2, line=2)
abline(h=0)
